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We analyze the phase diagram of the Bilayer graphene (BLG) at zero temperature and doping. 
Assuming that at the high energies the electronic system of BLG can be described within a weak 
coupling theory (consistent with the experimental evidence), we systematically study the evolution 
of the couplings with going from high to low energies. The divergences of the couplings at some 
energies indicates the tendency towards certain symmetry breakings. Carrying out this program, we 
found that the phase diagram is determined by microscopic couplings defined on the short distances 
(initial conditions). We explored all plausible space of these initial conditions and found that the 
three states have the largest phase volume of the initial couplings: nematic, antiferromagnetic and 
spin flux (a.k.a quantum spin Hall). In addition, ferroelectric and two superconducting phases and 
appear only near the very limits of the applicability of the weak coupling approach. 

The paper also contains the derivation and analysis of the renormalization group equations and the 
group theory classification of all the possible phases which might arise from the symmetry breakings 
of the lattice, spin rotation, and gauge symmetries of graphene. 
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I. INTRODUCTION 

Bilayer graphene^ (BLG) is a crystal which consists of 
two monolayers of honeycomb carbon lattice arranged ac- 
cording the Bernal stacking known from bulk graphite^l. 
In a Bernal stacked lattice, one out of the two sites on the 
upper monolayer resides directly over a site on the lower 
lattice, and the the other carbon atoms are on/under the 
centers of the hexagons (see Fig. Il| ) . Such a crystal has 
a very high symmetry with symmetry group 'D^d- 

This high symmetry may be lifted by the formation of 
correlated states of electrons. There is a plethora of ways 
the symmetry can be lifted, some of which have been 
discussed in the recent literature: the ferroelectric-layer 
asymmetric state^El^ the layer polarized antiferromag- 
netic stateSHB the quantum anomalous Hall state^ISElj the 
"spin flux"/ quantum spin Hall state^, the charge den- 
sity wave stat€p51 and an anisotropic nematic liquid^^^^. 
Some of the proposed phases above have a gap in the elec- 
tronic spectrum (ferroelectric, antiferromagnetic, spin- 
flux, CDW), whereas in the other phases (nematic, ferro- 
magnetic) no gap is formed. This large variety of possi- 
bilities makes the theory of electronic properties of BLG 
a very interesting and challenging subject. The com- 
plexity of the theoretical problem is compounded by two 
factors. One is a lack of precise information about the 
relevant interaction constants which determine the elec- 
tronic phase in undoped pristine BLG. The other issue is 
the competition between exchange energy contributions 
for a large number of candidate phases which makes the 
determination of the ground state non-trivial, even with 
precise knowledge of the interaction constants. 

On the experimental side, several contradicting obser- 
vations have been reported based on interpretations of 
the measured transport properties of suspended samples 



in terms of a gapful or gapless spectrum of electronic 
excitation4i2H16| However, all of these works as well as 
optical studies of BLGp21f22 indicate that the high-energy 
properties (but below Q.2eV) of BLG are well described 
by the two band model^ without interactions. This 
makes a comprehensive theoretical treatment of the prob- 
lem starting from the weak coupling even more timely. 
In this paper, we employ the previously developed RG 
approachPH to identify the possible scenarios of symme- 
try breaking phase transition in BLG at low temperature 
and zero carrier density. 

The tendency to form a state with spontaneously bro- 
ken symmetry is encoded in the system response to lo- 
cal symmetry breaking fluctuations, in particular in their 
mutual interaction. 
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Here 6py^ are operators creating local density fluctuations 
breaking lattice symmetry, with (5/5_4 = il;^ Mip expressed 
in terms of electron annihilation and creation operators 
ijj and ■0^, and gj^ are coupling constants. Each of the 
fluctuations Sp^ belong to one of the irreducible repre- 
sentations A (IrReps) of the symmetry group of the lat- 
tice. (A precise definition of the densities can be found 
in Sec.|H]). 

If Hamiltonian ( |1.1[ ) is dominated by one term with 
negative constant g^, we would expect it to be energet- 
ically favorable for a state with a non-zero expectation 
value of Spj\, to form, with the symmetry of the ground 
state determine by the corresponding IrRep, A. However 
if the coupling constant in the dominant term is positive, 
then the ground state is determined by the exchange en- 
ergy, which can be negative for not only for magnetic 
(ferro/antiferro) but also for non-magnetic orderings, be- 




FIG. 1. (Color online)Left panel: 3D view of bilayer 

graphene. The sites that sit on top of each other, connected 
by dotted lines, hybridize strongly and form bands with a gap 
of 7i ~ OAeV. The low energy electron live on the half of the 
carbon atoms that sit over/under the centers of the hexagons. 
Right panel: top-down view of the lattice. 



cause of the sublattice/valley matrix structure. Because 
of the large number of IrReps, this can result in a com- 
petition between many phases. Therefore, to determine 
the ground state of BLG we must know all the interaction 
constants g^ sufficiently well, especially when the dom- 
inant ones are positive. The situation is actually even 
more intriguing since attraction may result in a super- 
conducting phase with non-trivial Cooper pair structure. 
To add to the complexity of the problem, the values 
of the "constants" g^ axe not fixed. They change as 
a function of the energy scale £ within which the elec- 
trons establish the symmetry breaking correlations. The 
energy scale dependence, gyi(f ), may be calculated us- 
ing the renormalization group (RG) approach. In the 
RG approach the highest energy electron states are elim- 
inated and their effects incorporated into a redefinition 
of the parameters of the theory. The renormalization of 
BLG parameters starts at the energy scale 71 /2 w 0.2ey 
which limits the applicability of the two-band model 
with parabolic spectrum and initial conditions 3^(71/2). 
Then it is iterated until the lowest energy scale £ is 
reached. This energy scale £ is determined when the 
interaction energy in at least one of the the channels be- 
comes of the order of kinetic energy. After this scale is 
reached the mean field theory can be used to establish 
the electronic ground state. The necessary RG equations 
for the constants 5,4 and their interplay with Coulomb 
interaction. 
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have been derived for the full set of eight constants in 
Refs. illj (Similar in spirit treatment of Ref. 7 replaced 
Eq. ( |1.2| with the short range weak interaction.) 

Calculating g_4(7i/2) requires detailed knowledge of 
the microscopic orbitals which is not available at present. 
Therefore, in this paper we explore a wide variety of 
initial conditions (7^(71/2) for the RG to find possible 
electronic ground states for BLG. We can make some ar- 
guments to constrain the values of the (?_4(7i/2). The 
coupling gB2 which describes the interaction of dipoles 



oriented perpendicular to the bilayer (see Fig. l2| must be 
positive at high energy scales. The four "current-current" 
interactions g^^ ; 9Bi , 9Ei and g^" are only generated 
by virtual processes because of time reversal symmetry. 
Therefore we will set them to be zero at 71 /2. 

Also, it is interesting to note that in the value of 
5yl(7i/2) one has to take account of the interactions be- 
tween electrons via polarization of the lattice. Particu- 
larly, the in-plane TO-LO phonons at the F-point and 
TO phonons at the Brillouin zone corner have energies 
comparable to 71 /2, so that they mediate an attractive 
interaction via their virtual creation/absorption. These 
would give negative contributions to the bare values of 
gE2 and gE!-/ ■ Analogously, virtual LO—LA phonons from 
K - the Brillouin zone corners give negative contribution 
to the value of go- Therefore, we make no assumption 
about the sign of gE2 1 gs" 1 and go- A set of typical out- 
comes of the RG flow and the resulting electronic phases 
is shown in Fig. [31 

In Fig. [3] we reproduce the earlier reported resultfUn^l 
that for the initial choice oi g_A — the RG flow leads to 
a nematic phase. The nematic phase is a state with bro- 
ken rotational (but intact translational) symmetry cor- 
responding to representation E2 in Fig. (pi, mimicking 
the effect of anisotropic hopping along bonds with dif- 
ferent directions on the honeycomb lattice. This breaks 
the six-fold rotational symmetry by selecting one of axes 
of the lattice. In this state the electronic spectrum re- 
mains gapless but is significantly reconstructed from the 
unbroken symmetry state with two four-fold degenerate 
Dirac cones at low energy. The state has the same sym- 
metry and spectrum as uniaxial strain^'^, and we ex- 
pect that strain will, all else equal, favor the nematic 
phase. Figure [3] shows that the nematic phase is the pre- 
ferred ground state not only when 3^(71/2) — 0, but 
in a significant section of the g_4(7i/2) parameter space. 
In particular, the nematic phase always emerges from 
the part of the parameter space where bare electron- 
electron couplings causing intervalley scattering are zero 
{ga = gE'^ =9Ei' =0). 

In other parts of the parameter space explored in this 
work and illustrate in Fig. [3] the ground state appears 
to be anti-ferromagnetic (AF), with the Ai and B2 sub- 
lattices of two layers, see Fig. [l] are spin polarized in 
opposite directions. In the AF state the electronic exci- 
tations are gapped (though neutral spin wave excitations 
are gapless). Although the AF state prevails over a sig- 
nificant section of the parameter space, the combinations 
of high energy couplings which produce the AF state are 
not intuitive. For example, increasing the bare coupling 
gB2 does not necessarily introduce the AF phase. How- 
ever increasing the bare coupling gc makes the ground 
state AF. The reason for this counter-intuitive behavior is 
in the complexity of the RG flows. Since there are eight 
non-linearly coupled variables in the RG equation j^ l ^^ l, 
the RG flow is quite complicated, and the connection be- 
tween the couplings at low energy and the bare couplings 
at high energy is not obvious. 
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FIG. 2. (Color online) Sketches of the density and currents transforming according representations of the group T)'^^. In the 
case of a spin singlet symmetry breaking, the plus and minus signs represent charges, the blue lines persistent currents and the 
black bars represent bonds. The G, E2 and E'-l order parameters triple the unit cell; th e new Brav ais lat tice vectors are given 
by the dashed arrows. The representations are given in terms of Pauli matrices in Eqs. 12.111 and (2.101. 



Exploring a broader parameter space further we find 
more phases. A spin flux phase is found in a significant 
sector of the parameter space 3^(71/2), as seen in Fig. |3] 
This spin flux phase is a state with a persistent spin cur- 
rent circHng the honeycomb lattice rings, corresponding 
to the spin triplet form of representation Bi in Fig. |2] It 
may be viewed as the spontaneous formation of a strong 
spin-orbit coupling. It therefore leads to a gapped elec- 
tronic spectrum and possibly a quantum spin Hall effect. 

There are two more phases which appear to some de- 
gree in the phase space explored. One is a ferroelectric 
phase (FE). The FE phase a trivial band gap insulator 
where the bilayer becomes spontaneously charged like a 
capacitor. It is a completely gapped phase. It corre- 
sponds to representation B2 precisely the same represen- 
tation as AF but spin singlet, rather than spin triplet. 
Therefore, positive qb^ suppresses the ferroelectric phase, 
which appears in Fig. [3] only in the fine tuned corners 
corresponding to the applicability of the weak-coupling 
theory. 

We also found other phase a new superconducting 
phase (not shown in figure, see Fig. 9] for more details) 
which has the energy tantalizingly close to the nematic 
and ferro-electric states. It is a triplet superconduc- 
tor with a nontrivial Cooper pairing. Cooper pairs are 
formed between pairs of electrons with opposite valleys 
and opposite layer. The pairing is symmetric in exchange 
of valleys, but antisymmetric in exchange of layers. 

As usual, the singlet superconductivity appears only 
for the attractive interaction. From the first panel on 
Fig. [3] we see that it requires quite siginificant attraction 
in two channels gE"i9G < 0. 

Below we describe how the conclusions listed above 



have been reached. In Sec. |IT] we review the struc- 
ture of BLG, it's symmetry group and the low energy 
Hamiltonian. Section IIIII describes the resummation of 
the Coulomb interactions in the 1/N expansiorp32lH26]^ 
where TV = 4 is the degeneracy of the single particle 
spectrum. We then derive the RG equations that con- 
nect the couplings at low and high energy scales. In Sec. 
|IV| the results of the RG fiow equations are analyzed and 
augmented by a self-consistent mean field theory which 
produces a possible phase. Section [V| discusses the prop- 
erties of the emerging phases. In Appendix we describe 
the group-theoretic analysis of the phases of the BLG. 
diagram. 



II. MODEL 

The top view of the BLG lattice with Bernal stacking 
is shown on the right panel of Fig. [Tl Here we label the 
two layers 1 and 2 and the four inequivalent lattice sites 
Ai, El. A2, B2, with A2 directly over Bi. 

Calculation based on the minimal tight-binding model 
has established the following BLG band structure^. The 
A2 and Bi sites hybridize strongly and host states 
from the high energy bands with excitation energies 
> ji K OAeV. The low energy fermionic excitations in 
BLG belong to a four-component representation of the 
group 2?3d, exactly as in monolayer graphene. The four 
fermionic fields ip are conveniently joined into a 4-vector 
as follows 

v;*^{(^i^,^f)-^^(^f„-^ii,)^^}KK,, (2.1) 

where the ^s are true spinors including real electron spin. 
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FIG. 3. (Color online)Four cuts through the possible param- 
eter space of BLG. The predicted gap (or saddle point energy 
for only gapless nematic phase) is indicated by the color scale 
and the predicted phase is indicated. N is a nematic, AF is 
an antiferromagnetic phase, SF is a spin flux phase and FE 
is a ferroelectric. A flfth predicted superconducting phase is 
in a range parameters not shown, see Fig. [9] for more details. 
The g are coupling constants of BLG, with the subscript la- 
beling the irreducible representations in accordance with Fig. 
[2I and deflned in Sec. Ill] All boundaries are the first order 
phase transitions. 



This four dimensional space can be written as the di- 
rect product of the (AB) and (KK') spaces. We will 
use this to write all operators as the sum of direct prod- 
ucts T^^T^^ (Jc of Pauli matrices in each space. We 
define {t^^jT^^^ ,ai} as the Pauli matrices acting on 
layer, valley and spin, respectively, and define tq = 1, 
T± = {t^ ± iTy)/2. 

The symmetries of the BLG lattice consist of the two 
independent lattice translation ti and ^2, a C31 rotation 
by 27r/3 around one of the lattice sites; and two inde- 
pendent reflections: Rh, reflection across the y-axis, and 
Ry, reflection across the x-axis together with reflection 
through the plane midway between the graphene sheets, 
see right panel of Fig. [T] The reflections and rotations 
form the point group V^d- The groups V^d and Cqv are 
isomorphic and have precisely the same action on the 
plane. We also ignore the spin-orbit interaction which 
gives an additional SU{2) symmetry from the indepen- 
dent rotation of the spin. We will be concerned with the 



physics about K and K' points which are inequivalent in 
the Brillouin zone but are connected by Rh ■ Rather than 
dealing with two degenerate but inequivalent points we 
can triple the unit cell, which maps K and K' onto the F 
point. In this view, the point group D^d is expanded to 
■^3d ~ -^sd + ti^M + i2T^3d with the translation operator 
£1 with i\ = £2 and i\ = 1 (see e.g. Ref. E7)l . 

Ignoring the spin structure, the vector ^ transforms as 
follows under the action of the symmetry operators: 



ii^(r) . 


= e 3 ^. ^l; (tir) , 


Czi'ir) -- 


.e'rrr^^Car), 


R,Mr) -- 


-r^^r^^'i,{R^v), 


RvMr) - 


- T^'^rf^'^ (Ryr) , 



(2.2) 



There is also the time reversal synunetry operation given 

by 



V^^V^T; r^*f^^^ff^'a. 



A. Single particle spectrum 



We write the Hamiltonian for this model as 



(2.3) 



H = Ho + Hc 



H,,: 



(2.4) 



The single-particle part of the Hamiltonian in the two 
band modepJ reads (we will put ?i = 1 in all the subse- 
quent formulas) 



H. 
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KK' lAB ,2 
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.AB ,2 
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ij. (2.5) 



Here we ignore the "warping term"^ caused by the 
small skew hopping (73) since it would have a negligi- 
ble effect on the RG. We have deflned k± = k^ i: iky, 
and m = — 2"'''" 2 where 70 is the the interlayer inte- 

gral and tab is the interatomic distance. (Effect of the 
electron-electron interaction on the warping was studied 
in Ref. 11) The Hamiltonian in Eq. ( |2.5[ ) has the eigen- 
value spectrum 



e{k) = ± 
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2m' 



(2.6) 



where each branch is four-fold (spin and valley) degen- 
erate. The system described by Hamiltonian ( |2.5| has a 
higher symmetry than the underlying lattice. This larger 
symmetry is described by the SU{A) (E) U{1) group whose 
sixteen generators Mij are given by 



(2.7) 
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(2.8) 



An additional rotational U{1) symmetry extends the dis- 
crete rotation C3 to a continuous transformation given 
by V(^) -^ exp{-2iea^^)'ilj{R{e)r), where R{9) is the 
real space rotation by an angle 9 

The preceding discussion actually undercounts the 
symmetry algebra of the single particle Hamiltonian 
greatly, since they do not include the continuous parti- 
cle hole symmetry rotations^^. Including these rotations, 
the total symmetry group is Sp{8). However these extra 
rotations are not necessary for the following analysis. 
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B. Electron-electron interactions 



The Coulomb interaction 
He = - 



[(V'V)r (^^^)r,] 



(2.9) 



is the largest interaction energy in the system. The 
strength of Coulomb interaction on the length scale L 
is e^/L. The electron kinetic energy related to the 
same energy scale is l/(mL^) so the Coulomb interac- 
tion will dominate at the scale L = l/{me^), which is 
comparable to Bohr radius. However due to the gener- 
ation of electron-hole pairs, the Coulomb interaction is 
screened, leading to the reduction of interaction energy 
e^/L — )■ 1/(to7VL^). This screened interactions respects 
all the symmetries of the system and does not scale; 
therefore by itself it does not induce any spontaneous 
symmetry breaking of the lattice symmetry group. We 
will return to the quantitative description of the screened 
Coulomb interaction in Sec lHIl 

Any lattice symmetry breaking is captured by the scal- 
ing of the marginal short range interactions. These in- 
teractions also reduce the symmetries of the low energy 
model almost down to the crystal groupffiESIinEIl^ 

27r 
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(2.10) 



where have included a factor of — to make the couplings 
dimensionless. The V^^ symmetry of the two band BLG 
model forces various relations among the gij : 



9xx — Qxy — Qyx — 9yy — 9G 

9xz — 9yz = 9E2 ; 9zx — 9zy = 9e'^ 

9x0 = gyO = 5-Ei ; 9Qx — goy = gE{' 

9zo = 9bi ; goz = 9A2 ; gzz = gsa 



(2.11) 



Here we have labeled the couplings by the appropriate 
representation of ZJg^ schematically represented in Fig^ 
We note for future reference that the interaction terms 
gE2{'4'^'''xy^z^^ ^Y ^"^d gsii'ip^T^'^^p)'^ are invariant un- 
der the entire C/(4) (and which can be extended to 5*^(8) 
by including the particle-hole rotations^ symmetry of 
Hq. All other short range interactions, such as those of 
the form ~ [tp^^a^vil'vY oi' ~ IV'^V'JP can be always re- 
arranged into the form of H^^t by using standard Pauli 
matrix identity 2(5^j,^^'^' 
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FIG. 4. Definition of the elements of the diagrammatic ex- 
pansion. The thick line is the fermion propagator, the circle is 
the self energy from the single particle of the spectrum. The 
wavy line is the Coulomb propagator and the dotted line is 
the contact interaction. We separate the scalar contact inter- 
action goo from the other interactions. 



III. PERTURBATION THEORY AND RG 
EQUATION 

A. 1/N resummation 

For the Coulomb interaction we will use l/N as a small 
parameter, w here N = 4 is the number of degenerate 
fermion flavorgSEllESl, \\Ae achieve this expansion by per- 
forming the usual RPA resummation of diagrams (Figpl. 
Note that the coupling goo has the same matrix structure 
as the long range Coulomb interaction. We therefore re- 
sum the two together, i.e. we take the bare interaction in 
the RPA resummation to be 



V(")(g) 



27re2 



47r5oo 



Summing up the geometric series of terms in Fig. 
we arrive at the resummed propagator. 
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FIG. 5. Resummation of the strong Coulomb interaction in 
the 1/N approximation, a) Evaluation of the polarization 
loop, b) Definition of the resummed propagator, represented 
by the double wavy line. The scalar contact interaction is 
included in the resummation as it has the same matrix struc- 
ture as the Coulomb interaction, c) The non renormalization 
of the Coulomb vertex as a result of gauge invariance. {5Z is 
defined in Fig. Ima). 
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We further take the long wavelength limit, q ^ 0, 
where V^^^'lq)!! ^ 1. This gives us the approximate 
expression for the interaction propagator. 
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Since T) ex l/N we can use a perturbative expansion in 
1/A^. Note that we have neglected the higher-energy 
bands in considering the resummation of the Coulomb 
potential. However, the higher energy bands would only 
change the dielectric constant which cancels out of the 
final formula. 

Now we write the partition function as a path integral 
in imaginary time t over Grassman fields ij} and ij}\ 



Z = / D^jD^j^'e 



S 



j d\dtU^j^i>-H[ip\i;]\ 



(3.5) 



where H is defined in equation (2.4 1. Then, we perform 



the RG by integrating out all fermionic states with mo- 
menta q such that K > \q\ > K e~^ , where K is some 
ultraviolet cutoff regardless of lo. We will set Kq so that 
-?i'o/(2rn) — 7i/2, approximately the upper-limit of the 



FIG. 6. Diagrams included in RG equations, a) Single- 
particle particle weight renormalization; b) Renormalization 
of the mass m. c) Renormalization of the contact interactions. 



applicability of the two band model with the parabolic 
dispersion. We then rescale V' — ^ (1 + 5Z/2)'ip to keep 
the term /V'^;^'0 unchanged. This procedure has the 
benefit of not renormalizing the Coulomb vertex because 
of gauge invariance (see Fig. Isjc)). If we assign t an RG 
dimension 2 then at tree level the operator ip has RG 
dimension +1, and m, gij and the Coulomb interaction 
are marginal. 

There is a subtlety in the l/N treatment of the 
Coulomb interaction. Because of the behavior of the 
interaction in the limit g — > 0, w — >■ oo, some of the 
diagrams taken individually diverge faster than logarith- 
mically. For example, the self energy diagram (see Fig. 
loja)) gives the correction to the quasiparticle weight. 
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(3.6) 



where the variable x is defined by the substitution 
u! = xk'^/{2m). This integral is formally infinite since 
f{x) — > a; as a; — > oo. To understand this divergence, 
note that it conies from the region where the momentum 
k through the Coulomb line goes to zero. This corre- 
sponds to a spatially constant but time-varying poten- 
tial V{t). Such a potential is merely a constant shift 
in energy so that the Green's function are changed as 

G{ti,t2) ^ G (ti,t2) exp (ie f*^ V (t) dt\ . It is the sum- 
mation over the fluctuations of this phase that produces 
the divergence. However, in all observables gauge invari- 
ant quantities the fermion lines must come in closed loops 
which cancels out this phase, and so it cannot appear in 
any physical quantities. Reassuringly in all our calcula- 
tions this is the case. Indeed such a divergence cancels 
out from the correction to the electron mass. 
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Although the latter expression is divergent in the limit 
x ^^ cx). the sum, 



1 
2m 



1 



2m7V 



log K 



1^/(^)7^^, (3.9) 
27r [l + x^Y 



is convergent. Therefore the mass has a logarithmic de- 
pendence on cutoff, as expected. This enables us to write 
down the RG equations for the electron mass m. 



dlogm(£) 



'2N' 



£ = \ogiKo/K), 



(3.10) 



where 



ai 



1 

2^ 



2\3 



dxf{x){l-3x^)/{l + x^) 



-0.078. (3.11) 



Since ai/{2N) < 10 ^ is very small we shall neglect this 
mass renormalization for the rest of this analysis. 



B. Renormalization of the contact interactions 

We now consider the renormalization of the short- 
range interactions. Based on our assumption that the 
bare values g^ are small we will work to order g^ and to 
lowest order in 1/7V. 

The leading logarithmic corrections to the coupling 



constants of the contact interaction (2.101 are shown on 
the Fig. WicY Straightforward calculation of those dia- 
grams yielcPSl the set of RG equations for the 8 coupling 
constants g_4 

-df — ]^^(^2)., ^ g^, 

-i^-^B--2NA,gf^+±±C, 
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(3.7) where 
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(3.13) 
Here ^- is a sum over i,j = {0,a;,2/, z} excluding the 



combination i 



0,J 



and the summation conven- 



(3.8) tion is not used. The symbol 5{E2)ij is 1 when i = z 



and i = x,y and otherwise. By appearance there are 
16 equations contained in Eq. ( |3.12| . However several 
of these are identical due to the V^^ symmetry so there 
are only eight independent equations for the flow of the 
eight independent coupling constants. The numerical co- 
efficient ai is defined in Eq. (|3.11 1 and 



Ck2 



aa 



dx 2f{x) 
2^(14-2:2)2 

dx fjxf 
27r 4(1 + 2:2)2 



.469, 



.066. 



(3.14) 
(3.15) 



The term 2NAijgfj in Eq. ^.12| corresponding to 



leading loop diagram (v) in Fig. fol cj is naively the most 
significant quadratic term in Eq. (3.121, because it is 
leading in A^. This term represents screening of repul- 
sive interactions in the charge channel as expected in a 
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fermionic system (since Aij > 0). Note that this term 
is actually zero for the representation £'2 and A2, be- 
cause these interactions commute with the single par- 
ticle Hamiltonian. Therefore, to lowest order in 1/iV, 
the interactions E2 and A2 are unscreened and free to 
grow strongly attractive. We hasten to add that that 
the higher order terms in Eq. (3.121 are very important 



and one cannot understand the behavior of the RG flows 
based only on the leading terms. 

The single particle Hamiltonian is off-diagonal so there 
is a contribution of order g^N~'^ to the coupling f;^^ from 
the two Coulomb line diagram in Fig. l6jc)(iii,iv). These 
may be calculated, 
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(3.16) 
From this, it follows the that the free field point gA — ^ 
is not a fixed point. Even if the the system starts with 
all bare couplings (7^(71/2) = it will fiow under RG to 
have finite gE2 and gB^ with the other couplings fixed to 
zero by the SU{A) symmetry of the single particle Hamil- 
tonian. To demonstrate the behavior in this regime we 
ignore momentarily g^i which gives us a single equation 
for gs^, 
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(3.17) 
Since the first term on the RHS is negative, there can be 
no fixed point and gE2 flows to —00 regardless of the ini- 
tial conditions. (This holds whether we treat Eq. (3.17) 
to lowest order in N or simply plug in A^ = 4). Accord- 
ing to the mean-field theory (see Sec 
a nematic ground stated with transition at 



IV I, this suggests 
lOOmK. 



C. Applicability of our approximations 

Let us turn to the justification of only including the 
diagrams Fig IgJc) in our treatment. Notice that it is 
different from the conventional 1/A^ approximation, see 
e.g Ref. 1311 There are two issues: (1) there are two loop 
diagrams which are leading order in 1/N but are not 



included. Fig. ItJc); and (2) there are diagrams that are 
subleading in N which are taken into account (compare 
the bubble and ladder diagrams in Fig. l6jc). 

To address the first issue, let us discuss diagrams of 
Fig. [Tie in more detail. They are non- vanishing only for 
gE2 and have the form SgE2 ^ ^j^ (Cln [K) + Z?ln(fc)). 
The term In is produced by two iterations of the the RG 
equations (by substituting diagrams (iii) and (iv) into di- 
agram (v)) but the second term does contribute to the 
linear term in the RG equation 2a2AE2 -^ 2a2^B2 + ^■ 
The constant D, however, depends on the cutoff scheme 
so that the term linear in gE2 in the RG equation for gE2 
is not known (for the other constants it is well defined). 
Fortunately, is does not matter for the divergent behav- 
ior at large N. Consider the situation with all other con- 
stants except gE2 fixed to zero, keeping only coefficients 



leading in 1/A^, compare Eq. (3.171 



dg 



E2 



as 



a\ -I- 2a2 

aT 



-9E2 



2iVg|^ 



(3.18) 



The quadratic term dominates the constant term when 
9E2 ^ N~^ at this point, but then the linear term is 
smaller by a factor of \/^/N <C 1. Thus, contrary to 
initial appearance, the linear term is of higher order in 
1/N for gE2 - so that we leave it in Eq. (3.121 only 
for simplicity. It makes essentially no difference to the 
evolution of the RG equations. 

To address the second question we notice that the bub- 
ble diagram Fig. [6fc(v) contains an extra factor of N in 
comparison with diagrams (vi, vii,viii). The latter dia- 
grams are not diagonal in terms of the coupling constant, 
as given by the tensor C^'5„„, whereas the bubble dia- 
gram is oc Ngf;: by construction. The large amounts of 
constants involved in the non-diagonal term may over- 
come the factor of N in the diagonal terms; therefore 
keeping both is legitimate. The higher order terms may 
be considered as 1/A^ corrections to the tensors Aij and 
^kimn respectively. For example. Fig. ]J\&) is a leading 



t 



1/N correction to Aij, whereas Fig. ujh) is a leading 
1/N correction to C^^;^„, even though the two diagrams 
do not have the same in order in N. 

Finally, we compare our treatment to the existing the- 
oretical contributions. The work of Vafek and Yangl^l 
is similar in spirit but contains only the Gi and B2 
out of the eight possible representations and treats the 
Coulomb interaction as short range. The later work of 
Vafek^contains the RG equations for the full eight con- 
stants but again treats the Coulomb interaction as short 
ranged, and does not attempt to describe the general 
structure of the RG flow. The treatment of Ref. ^ is 
completely at the mean-fleld level and corresponds to 
counting only the diagrams from Fig. pf marked (i) and 
(ii), which is not a parametrically justifled approximation 
as well as considering only the B2 representation. The 
RG equation of Ref. E] essentially ignore the spin- valley 
structure but still does not retain the necessary number 
of coupling constants to describe even that limited situ- 
ation. 
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FIG. 7. a) Schematic representation of the "bubble" diagrams 
where the shaded blob represents all possible connected dia- 
grams, and arbitrary Coulomb propagators may be added, b) 
Similar representation of the "ladder" diagrams. The leading 
diagrams from both of these groups are included, even though 
this is not strictly parametrically correct, c) A second loop 
contribution to the anomalous dimension of the coupling gE2 
which is disregarded. 
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FIG. 8. (Color online)Plots of the coupling constant as a 
function the running RG scale £ — \og{Ko/K). There are 
eight running couplings labeled by the corresponding repre- 
sentation. The density-density couplings are given by solid 
lines, the current-current couplings by dashed. The graphs 
end when the couplings become of order 1/N = 1/4. They 
reach a singularity a finite £ soon after the graph ends. 




IV. RG FLOWS, THEIR TERMINATION AND 

RENORMALIZED MEAN FIELD TREATMENT 

OF SYMMETRY BREAKING 



In this section we describe the numerical analysis of the 
coupling constant RG flows described by Eq. (3.12) and 



show that there are no weak coupling fixed points. The 
divergence of coupling constants in 2D at zero tempera- 
ture indicates spontaneous symmetry breaking. (Unlike 
in ID the quantum fluctuations in 2D are not infrared di- 
vergent and do not destroy zero temperature phases.) We 
analyze the resulting phases within mean-fleld theory, us- 
ing the coupling constants renormalized by the RG. This 
treatment is superior to simply doing mean fleld start 
from the high energy scale since in that case the large 
logarithms are not summed in a controlled fashion. 



A. General structure 

If the initial RG conditions are such that 5^(71/2) ^ 0, 
then the 5[/(4) symmetry is absent and we must consider 
the flow of all the coupling under the RG. Determining 
if there exist any fixed points cannot be done analyti- 
cally as it requires solving a polynomial of the 64th order. 
However, a numerical solution shows that there exist no 
fixed points. Therefore at least some of the couplings 
must grow infinitely. At the same t ime, t he leading term 
quadratic in the couplings in Eq. (3.121 ^ Ng'^, always 



with a non-positive coefficient, so we expect generically 
that large positive coupling to be driven back to zero. For 
all positive initial coupling constants, this means that the 
system will be driven to the free field point until g^^ be- 
comes large and negative. This behavior is confirmed by 
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the numerical evolution of the RG equations (see Fig. [s] 
where gE2 always becomes negative and increases until 
the other couplings diverge). 

Note that although we have set the current-current 
couplings gBi , 5^2 1 9Ei and gE" to zero initially, they are 
generated through renormalization. The examples of the 
RG equation shown in Fig. [8] indicate the current-current 
couplings become of the same order as the density-density 
couplings at low energies. Therefore we cannot ignore the 
current-current interactions when analyzing the ground 
state, and ignoring them would lead to misleading results. 

The coupling g^^ lias been given special emphasis in 
some of the earlier studies^. We find that in the RG equa- 
tions it does not seem to play an exclusive role, as can 
be seen in Fig. (l9|, where it is screened efficiently - the 
leading term in the RG flow is 2Ng^ . Note that large 
positive initial gs^ does not provoke a phase transition to 
the AF state on its own (see Fig (l8|). Once g^^ becomes 
relatively large the presence of a finite g^a will change the 
structure of the fiow, especially since it breaks the 5C/(4), 
however not in a marked way. For example, starting with 
all other couplings set to zero except for y^^, gE2 still 
becomes the most significant negative coupling, and the 
nematic phase is the preferred phase. As a result, g^^ is 
perhaps the least important of the four couplings. This is 
not a conclusion that can be reached on general grounds, 
but only by solving the detailed RG equations over a 
broad range of parameters. Moreover, at least some of 
the couplings behave non-monotonically. Initially negli- 
gible coefficients may end up diverging quickly (e.g. 17^12 
in FigllsFa)). At the same time a couplings that is not 
large at the end of the RG flow may change the character 
of the flow in the initial stages. 

The RG equations contain terms up to second order 
in g. Therefore, it may be easily seen that since there 
is no fixed point, the couplings always go to infinity as 
9 A ^ ^a{£o — i)~^ where £q gives the value of the singu- 
larity in the RG flow and the A_4 determine how quickly 
each constant diverges. Mathematically, there are six sets 
of {A^} that satisfy the RG equations and are stable to 
perturbation. One might be tempted to determine the 
ground state, using this mathematical feature, via the 
coefficients A^. However these coefficients are meaning- 
ful only for the RG at g_4 ^ 1 which is outside the range 
of validity of the proposed theory and the RG equation 
(3.121. Moreover, g — >■ 00, indicates an instability to- 



wards a broken symmetry states, so that we shall use a 
mean field theory starting from the energy scale where 
some of the couplings become sufficiently large. 



B. Ground state energies within renormalized 
mean-field approach 

The unbounded growth of coupling constants in the 
RG flow generally indicates the development of a spon- 
taneous symmetry breaking and the opening of a gap. To 
describe the corresponding phase transition we use a self- 



consistent mean-field theory. The self-consistent mean- 
field theory is implemented by replacing all possible pairs 
of fermions in the quartic interaction terms with their 
mean values. For this we introduce the Gorkov-Nambu 
vectors which adjoin the two 8-component vectors ip and 
ijj^ as follows. 



*wHr(^U))' 



N 



(4.1) 



N 



where T = ir^^ '''y^'^y is the time reversal matrix. 
We also introduce Pauli matrices t^^ ^ acting on the 
Nambu space. The vectors "^^ and "^ satisfy the condi- 
tion 

¥{k) = i-^\-k)T^f (4.2) 

We rewrite the interaction terms with the help of the 
Nambu vectors 



^(V.Vf^"r/^V^)^ = jE.9.(*t.M..* 



(4.3) 



Here KL 



-KK'AB^N, 



r"" Tj T^ ai acts on the 16 dimensional 
space spanned by the Nambu vectors and we write s as a 
shorthand for the list (ijkl) . The couplings gs are defined 

as gabzQ = gab, ffoozO = 9oo, 5o000 = 9a0 and goboo — gob, 
where a,b — x,y, z and all other constants are zero. The 
factors of t^ are necessary since we must have 



i^r^'T)Ml{^r^'T) 



-Af, 



(4.4) 



to satisfy both Eq. (4.2 I and fermion anticommutivity. 



The mean field approximation consists of replacing 
pairs of fermionic operators in Eq. (4.3 1 with their ex- 
pectation values as follows: 



s s \ 



^T.M'^.^^T.^f's.^ 



-2*^- (Ms- (*«)*^)-M,') *- - d^-M" 



■^^y 



tr 



((*(8)*^) 



•M. 



(4.5) 



Here we have used the fact that, according to Eqs. (4.2 1 



and (4.4 1 



to combine the Cooper and Fock terms. 

Now we assume that the there is some nonzero expec- 
tation value of the fields which corresponds to a non-zero 
order in one of the phases classified in the Appendix. 

i^^^')^-^^Oi:^n^% (4-6) 
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FIG. 9. (Color online)Plot of the broken symmetry phase and energy scale for BLG as a function of bare coupling constants. 
When the energy scale is less EhiTr ~ ImeV the broken symmetry state is in competition with the Lifshitz transition and the 
BLG may remain metallic (with eight Dirac points). Dashed line indicates the region where the triplet superconducting phase 
(SC) is very close in energy (but slightly above) to normal (nematic or ferroelectric) states. 



where matrices Af^ are specified for each phase A in the 
Appendix and A^ = {A^} is the order parameter, which 
is singlet or multi-component depending on the phase. 
Below we will use the notation |A_4p = J^a I^SP- 
The effective interaction constant c^ are defined for each 
phase as, 

1 



CA 
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9s 



SsA - 



4iV2 



tr 



Ms Ma 



(4.7) 



The assumption of a finite expectation value is consistent 
only if CA < 0. The interaction mean field energy is 
therefore. 



^--2E^S(*^-^^5-* 



87rc_4 



(4.8) 



Including the effect of the single-particle Hamiltonian Hq, 
defined in Eq. (|2.5|, the total mean field Hamiltonian is. 



HMF = lY.'^Hk) 



T^orf+^A^M^ 



*(fc) 



Sttca 
(4.9) 



for fixed values of the order parameters A|^. 

In the spirit of the Hartree-Fock or BCS theory, we 



diagonalize Eq. (4.9 I to obtain the ground state energy 
per unit area, 

Emf{{Aa}) = -N f ^eik,AA)- 

J^<£, 27r 



87rc^(£c) 
(4.10) 
Here e(fc, {A_4})) are the positive eigenvalues of the ma- 
trix Hqt^ -\-'Y^^/^'ji^M'^ and the factor of N comes from 
the degeneracy. The energy scale £c is the energy scale 
at which we stop the RG. The integral in Eq. (4.101 
evaluates to 



const 



-<£c 



27r 



a A + /3.A In ■ 



Ml 



(4.11) 
where a^, and /3_4 are coefficients that depend on the 
phase A. The coefficients a a may be explicitly calculated 
from Eq. (4.101. We list here a a for the states where 
A • M may be written as X^jj/t; ^i'^^j'^kXir^^Tj^^ ^kT^^ ■ 
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Because of the symmetry of Hq there are only three in- 
dependent coefficients. Labeling the coefficient a by the 
representation and using the superscript n or s to denote 
normal or superconducting we have 



»M = «A. 



ar 



ar 



a. 



= 1 

1 

^2 



-+log2w 1.19, 



a 



ar 



"G 



ar 



E, -C.S, -aG= 4 +log 2-0.94. 

(4.12) 
The coefficients for the spin singlet and spin triplet nor- 
mal states are the same because of symmetry. These co- 
efficients are sufficiently close to one that w e have simply 

taken a^ ~ 1. the term /J^ln^^ in Eq. (4.111 should 
be interpreted as the continuation of the KG how from 
the scale £c down to the energy |A|. Although obtained 
by using mean-field theory, since the flow of c^ is gov- 
erned by the RG equation (3.121, we have to replace the 
logarithmic correction to Eq. (4.11 1 with the evaluation 



of cj( at the energy A. The mean field energy density is 
therefore written as. 



E{Aa) = 



Stt 



-2N - 



1 



Ml 



c^<0 



ca{\^a\)_ 

(4.13) 
The ground state may now be determined by mini- 
mizing Eq. (4.131 with respect to A_4 with the c_4(A) 



obtained by numerical integration of the RG equations. 
The ground state energy gap will then be equal to the 
value of |A| at the minimum. 

It is important to note that we expect to find this 
minimum when the coefficient c(|A|) ^ -^ which is in- 
side of range of validity for our RG equation, g — 1/A^. 
The remaining subtlety is the inclusion of the long- 
range Coulomb interaction into the mean-field descrip- 
tion. Usually, it enters in the statically screened limit 
(7oozO — 1/-/V and does not diverge at the transition. 
Furthermore, according to Eq. (4.7), this constant can 



produce only finite l/N"^ correction (positive to all su- 
percoducting states and negative for all normal states). 
Therefore, we will neglect goozo in the further manipula- 
tions. 



THE PHASE DIAGRAM 



^LiTr ~ ImeV. In this case the renormalization of cou- 
pling constants is stopped at EliTti spontaneous symme- 
try breaking does not occur and the system remains in 
the symmetric state with the four Dirac cone spectrum. 

Figure [9] shows the results of the RG analysis in 
terms of the resulting symmetry broken phases. We 
find that there is significant variation in the scale £, 
\og{£Q/£) '^ 1 -^ 20, as expected from the wide range of 
couplings analyzed. If we consider small initial couplings 
the (7_4 <g 1/N then the RG is driven by the constant 
term and log(i5o/^) — 10 irrespective of the initial condi- 
tions, resulting in a symmetry breaking only at extremely 
small energy scale £ ~ 10~^meV. 

There is a variety of phases that have been proposed as 
the ground state of BLG that we do not find. An anoma- 
lous quantum Hall state (QAH) state was suggested in 
Ref. Ss corresponding to the representation g^^ with or- 
der parameter {ip'r^ tjj). This is not found as a ground 
state in our analysis. In the same paper, and in Refs. [5] 
and 1321 a large manifold of quantum hall ferromagnetic 
states were suggested containing the representations £'2 , 
A2 and B2 in both spin singlet and spin triplet represen- 
tation. All of these states were considered as degenerate 
appealing to the S'C/(4) symmetry of the single particle 
Hamiltonian, Eq. ( |2.5| . In both cases the artificial S't/(4) 
symmetry was assumed to be exact , which is contradicted 
by the importance of the short range interactions we find 
here in solving the complete set of RG equations. 

In subsections below, we discuss the details of each 
phase. We will present a comparative flows of the cou- 
plings defined in Eq. (|4.7| to illustrate the competition 



between phases, see Figs.|10|-[T5 



A. Nematic Phase (N) 

In the nematic phase, there is a finite expectation value 
for the order parameter {^^t^-^ T^y^), breaking the ro- 
tational symmetry of the system from six-fold to two-fold 
while maintaining translational symmetry. The order pa- 
rameter is characterized by enhanced electron hopping in 
one direction. The interaction energy for this phase de- 
pends on the combination of parameters obtained from 



Eq. (4.71 



1 



1 111,,,, 

j9e^ +9e2 + ^5Si - ^9e'{ - ^9A2- (5.1) 



The result of minimizing Eq. (4.131 is presented in 
Fig.[9j We find by extensive numerical investigation only 
five out of the possible sixty-four phases enumerated in 
Appendix (10 in the charge channel, 22 in the spin and 
32 in the Cooper channels). They are the nematic phase, 
the antiferromagnetic phase, the spin flux phase, and in 
the corners of the parameter space of bare interaction, 
ferroelectric phase, and singlet and triplet superconduc- 
tor phases. 

Note that it is also possible that the resulting gaps 
are smaller than the energy of the Lifshitz transition^^. 



It is the preferred ground state in the absence of in- 
tervalley scattering, and it is also generally the ground 
state when the bare gE2 coupling is negative. Note that 
a negative contribution towards bare g^^ comes from the 
electron-electron interaction via the polarization of the 
lattice, ie. via virtual excitation/absorption of in plane 
phonons near the F point. The nematic phase is also 
the ground state over other large parts of the parameter 
space as can be seen in Fig. [9] This reflects the fact that 
the coupling gE2 almost always becomes negative rapidly 
(see Fig. [§. 
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FIG. 10. (Color online) Interaction energy as a function of 
energy scale for selected phase. For these initial conditions 
the nematic phase is the ground state 



FIG. 11. (Color online)The interaction energy as a function 
of energy scale for some choice of initial parameters. In this 
case the ground state is the AF phase. 



We previously proposed the nematic phase as a pos- 
sible ground state in based on a more limited analysis 
of the RG equation&ll- Using a similar renormalization 
group analysis Vafek and Yang^ ^^ also find the nematic 
phase as a possible ground state, supporting the analysis 
in this paper. 

The most notable characteristic of the nematic phase 
is that it remains gapless, but with the parabolic bands 
reconstructed into two Dirac mini-cones at energies less 
than I A|. The nematic phase would show metallic behav- 
ior in conductance measurements, but decreasing density 
of states at low densities. The nematic state preserves 
time reversal symmetry and so it should not show asym- 
metry between positive and negative magnetic fields in 
magneto-transport. The nematic order parameter trans- 
forms in the same representation of D^d" as uniaxial 
strain so that strain will couple directly to this state. 
It is possible that strain could induce a transition into 
the nematic phase^, even if unperturbed BLG chooses 
another phase as the ground state. 



B. Anti-ferromagnetic Phase (AF) 



The anti-ferromagnetic phase is defined by a nonzero 
expectation value of {ip^T^^r^^ <^^j}- This state cor- 
responds to opposing magnetic moments on the A and 
B sublattices. The orbital part breaks the refiection 
symmetry between the two sublattices but otherwise pre- 
serves the I?3^ symmetry of the BLG in its entirety. 

The exchange energy depends on the following combi- 
nation of coupling constants obtained from Eq. (4.7 1: 
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:9G 



+ 4 i9E^ 






X 1 

9Ei +gE['j - gffisi - 



;9A2 



(5.2) 



The AF promoted strongly by the coupling ga, with the 
factor of four coming from the dimension of the repre- 
sentation G. This effect is amplified by the sensitivity of 



the RG equations to the coupling gQ. As a result even 
small values of go near the free field point make the AF 
state the ground state. The AF state is also promoted by 
negative gE" and is generally the ground state when we 
start with negative gs"- Again we emphasize that these 
conclusions come from the combination of RG equations 
and interaction energy, not just from the interaction en- 
ergy alone. Other authors have proposed this AF state as 
a possible ground state. Kharitonov- suggested the AF 
state based on experimental evidence and simple mean 
field theory arguments applied at the high-energy scale 
directly. The structure of our RG equations indicate that 
such a simple mean field theory is not applicable to BLG, 
although it does suggest the same state. 

The AF phase is expected to be an insulating state 
with activated gap behavior in transport measurements. 
Although the magnetic field does not couple directly to 
the AF state since it is antiferromagnetic, the Zeeman 
energy splitting does break the SU{2) spin symmetry 
of the system (spin-orbit is negligible). Therefore the 
gap should show a decline in tilted magnetic field. The 
LAF state is adiabatically connected to a quantum Hall 
ferromagnetic state at higher magnetic field. A lack of 
features between the zero and high magnetic field state 
might be the evidence that the zero magnetic field state 
is AF^. 



C. Spin flux (spin Hall) Phase (SF) 

The spin fiux (SF) phase (elsewhere called a quantum 
spin Hall state) is defined by the finite expectation value 
of (■(/; V^^CTi/;) . The effect of this on the electrons is equiv- 
alent to the development of a finite spin orbit coupling, 
and gaps the electronic spectrum. It may be viewed as a 
state where spin currents circle the honeycomb rings, or 
as a quantum anomalous Hall effect state, but with op- 
posite signs for opposite spins, producing no net charge 
current, so that this state preserves time reversal invari- 
ance. 

The interaction energy of the SF depends on the com- 
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FIG. 12. (Color online) Interaction energy for selected phases 
a function of energy scale. In this case the ground state is the 
spin flux state shown in orange 



bination of coupling constants obtained from Eq. (4.7) 



1 



1 



CSF 



{gs'^' - gE2 + 9Ei - gE[')~^ {9B2 + .gsi + .9^2) • 

(5.3) 
The SF state has a similar matrix structure as the AF 
state and therefore a similar interaction energy. However 
it is not promoted by large go, unlike the AF phase. 
Therefore, geirerically, large ga generally suppresses the 
SF in favor of other states as can be seen in Fig. [9) 

Analogously to the case of spin-orbit coupling in mono- 
layer graphene^^ , the finite value of the spin flux OP may 
create a "spiir Hall effect" with quantized spin Hall coir- 
ducivity The edge states and insulating bulk imply a 
quantized conductance of 4e^//i, unless they are local- 
ized by magnetic and intervalley-scattering disorder. The 
state is time reversal invariant, so no transverse conduc- 
tance at zero magnetic fleld is possible. 



D. Limits of applicability: Ferroelectric Phase 



FIG. 13. (Color online)The evolution of the interaction 
energy of the selected phases as a function of the energy scale. 
In this case the ground state is a ferroelectric phase. The 
difference in energy between the F phase and nearby phases 
is quite small. 



The FE state was also proposed in Ref. ^ However 
this analysis was based on a flawed mean field theory 
treatment counting only the diagrams from Fig. m\ (c) 
marked (i) and (ii) as well as considering only the single 
parameter scaling theory with one interaction constant 
in the B2 channel. These are not parametrically justified 
approximations . 

The ferroelectric phase is a completely gapped state, 
with neither neutral nor charged excitations. It is also a 
trivial insulator in that it does not possess protected edge 
modes. Therefore, it should display insulating transport 
behavior. An external field perpendicular to the BLG 
flake would promote the FE phase, increasing the gap 
due to the interlayer symmetry breaking at the single 
particle level^. This does not s eem to take place in any 
of the recent experimentaSHSl on BLG, where the ex- 
ternal transverse field destroyed the zero-field state, and 
introduced a distinct state determined by the interlayer 
asymmetry. 



The ferroelectric (FE) phase is characterized by a non- 
zero expectation value of {ip^Tf^T^^ tp). It is a sponta- 
neous charging of the BLG with opposite charge on the 
two layers. 

The interaction energy for this phase depends on the 
combination of couplings obtained from Eq. (4.7 1: 
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^ [9e'^' ~ 9E2 + 9Ei - 9e[') - g5Bi - g5A2 ■ 

(5.4) 

We find that the ferroelectric phase is strongly suppressed 
by positive bare ^s^ and the development of the FE phase 
requires the ^s^ coupling to diverge to negative infinity, 
which can only happen in a small sliver of the phase space 
where the combination of the higher order diagrams con- 
spire to drive gB2 negative. Even in this section the en- 
ergy difference between the ferroelectric and competing 
phases is never large, and it may be that in a more ac- 
curate calculation it never appears as the ground state. 



E. Limits of applicability: Superconducting Phases 

The singlet superconducting (SS) phase is character- 
ized by the usual order parameter (■0^T^^). The cou- 
pling for this phase is found from Eq. (4.7 1 to be: 



CSS 



1 



9G 



-,9B2 



+ T y9E2 + 9e'^ - 5Bi - 5bi', 



-,9A2 



>9Bi- 



(5.5) 



The stability of this phase requires very significant nega- 
tive couplings from the very beginning and the phase can 
not arise from the purely repulsive interaction. 

The triplet superconducting (SC) phase has the order 
parameter {tp^T^^ al'if/^), with the pairing function of 
the opposite signs in the K and K' valleys. 

Its interaction depends on the combination of couplings 
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FIG. 14. (Color online) Interaction energy for selected phases 
as a function of the energy scale. In this case the singlet 
superconducting phase is the ground phase. 




^ E(meV) 



FIG. 15. (Color online) Interaction energy for selected phases 
as a function of the energy scale for initial repulsive short 
range interaction. In this case the triplet superconducting 
phase is very close to the ground phase. Even though csc 
appears to be the most negative, minimization of Eq. (4.131 
gives the nematic state as a preferd ground state. 



considering the high energy couplings in detail, because 
different couplings lead to different ground states. More- 
over naive expectations about the importance of certain 
couplings are not borne out and all plausible combina- 
tions must be considered. In particular excluding inter- 
valley scattering leads to misleading results. The previ- 
ously reported nematic state is the ground state for a sig- 
nificant fraction of these couplings. Of the large number 
of ground states that are possible we find only five ap- 
pear. They are nematic, antiferromagnetic, ferroelectric 
and triplet superconductor, as well as a "spin flux" phase 
not previously proposed. The nematic, antiferromagnetic 
and spin flux phases seem the most likely candidates. 

The present work may be extended in a variety of ways. 
The accuracy of the renormalization group equations may 
be improved and validated by considering higher order 
diagrams and a more detailed mean field theory con- 
structed. One may also try to connect the value of the 
couplings at the scale 71 /2 sa 0.2meV with their value 
at the bandwidth of the tt orbitals. The possible phase 
transitions between the proposed states and the behavior 
of domain walls between regions of different phases may 
be needed to properly account for the transport data. 

Such improvements aside, the unique challenge of theo- 
retically determining the electronic ground state of BLG 
has been laid out. The problem naturally involves the 
competition of an uncommonly large set of phases and 
interactions. Truncating the theory to a more tractable 
subset does not appear to give accurate results. Instead 
the problem must be attacked in its full complexity. 
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obtained from Eq. (4.7 1: 
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For the repulsive interaaction, the triplet SC phase only 
appears as a stable phase at large couplings. 

To conclude, both superconducting phases appear only 
on the limits of the applicablity of the theory. More- 
over, finite value of the underscreened Coulomb interac- 
tion push the energies of those stated further up. It is 
therefore possible that the appearance of the SC phases is 
merely an artifact of our inability to deal properly with 
strong couplings, and that the superconductivity does 
not belong to the actual phase diagram. 



VI. CONCLUSION 

Several results have been established in this paper. 
The ground state of BLG cannot be understood without 
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Appendix: Group theory for phases in BLG 

In this appendix we classify the possible phases of 
BLG. We will use the matrix notation defined in Sec. |TT1 

Phases are defined by all possible expectations A = 
(^ (g) ^^) that belong to an irreducible representations 
(irrep) of the symmetry group Q of BLG. Every phase de- 
fines a subgroup % oi Q consisting of all operations that 
leave A invariant. Two phases within an irrep are distinct 
if their invariant subgroups are not conjugate. (Recall 
two subgroups % and %' of a group Q are conjugate if 
they is an element oi g & Q such that g'Hg~^ — T-L'). This 
definition is correct is in the sense that it gives all physi- 
cally distinct states that may be reached via a second or- 
der phase transition at the highest critical temperature. 
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per the usual Landau theory. 

Let us notice, however, that the anomalous averages 
belonging to the same irrep of the original group G 
may correspond to the different phases. For example. 
A, = t^^'t^'"' and A„ = t^^t^""' + rf^r^''' are 
both charge density waves with a tripled unit cell trans- 
forming in the G representation. However, these two 
phases are distinct since A// is invariant under rotations 
by 27r/3 around a lattice site, whereas A/ is not invariant 
under any conjugate operation, see Fig. [16] As a fur- 
ther example, a canted anti-ferromagnetic phase would 
be given by A = t^^t^^ Ox^'^z- This is not considered 
in the present classification since it does not belong to an 
irrep of Q (it is a linear combination elements of the B2 
and Ax representations). Of course all such mixed states 
may be constructed from linear combinations of phases 
in our classification. 

We classify the phases according to the symmetry 
group Q of BLG which is effective at intermediate en- 
ergies, see Sec. [IT] In this regime the effect of RG irrel- 
evant perturbations such as Umklapp scattering may be 
ignored, and RG relevant but weak perturbations, such 
as trigonal warping and spin-orbit coupling, may be ne- 
glected. This approximation should be effective at ener- 
gies between 71/2 = QfleV and Elitt = ImeV, which 
contains any energy scale associated with spontaneous 
symmetry breaking. In this regime the symmetry group 
is 

I?lrf ^ X I?f*r"^ X 5C/(2)("P™) X [/(l)(9™9'=) X T. 



inf 



inf 



(A.l) 

is generated by infinitesimal 
spatial rotations C and inversion Re of the BLG plane. 
The second subgroup, 2?i,i™" , is generated by an in- 
finitesimal translation t and refiection Rt- The action 
of these operators on the low-energy electrons is given in 
terms of the Pauli matrices: 



The first subgroup, "D^^f 



C 



^trf^ 



Rt 



Rc ^r^^rf^' 
Rt ^r^r^' 

. n AB KK' 

±ty — I y I y 



(A.2) 



Note that these two groups commute with each other, 
unlike true translation and refiection. In the presence 
of the appropriate symmetry breaking these are reduced 
down to the D^d" group discussed in the text. In this case 
the continuous translations and rotations exp {9tt + 0^0) 
become discrete with Ot^c = 0,±27r/3 and the inversions 
Rt and Rc become Rt ■ Rv and i?„ of Sec. [ll| respectively. 
The S'C/(2)(*P™^ is the group of spin rotations, which is 
decoupled from the physical rotations because there are 
no spin-orbit interactions. It is generated by rotations 
S = (Sx, Sy, Sz)- The high symmetry axis of any phase 
will always chosen to be z. Infinitesimal rotations around 



the z axis are represented by Sz- There are also refiec- 
tions and inversion of the spin space, but these do not 
distinguish any phases, so we will suppress them. 

The gauge groups acts by multiplication. We label the 
infinitesimal generator g, which acts on the wavefunc- 
tions ijj simply by gip — i-tp. 

The time reversal operator T commutes with all of the 
above except the gauge generator TgT — —g. It is given 
in terms of Pauli matrices by 



T; 



KK' AB T^ 



(A.3) 



where K is complex conjugation. 

This defines the symmetry group Q fully. A further 
enlarged symmetry group isomorphic to f/(l) x f7(4) is 
considered at points in the body of the text [see Eq. ( [2.7[ )[ 
and other works, but there is no reason to expect that this 
will every be an accurate approximation. We will proceed 
to categorize the phases according to the symmetry group 
Q. One may always collapse classification to obtain the 
distinct phases under the artificially enlarged symmetry 
groups. 

In the following subscections we present the tables enu- 
merating the possible symmetry breaking for the singlet, 
triplet, and superconducting phases each. Each table 
is split into subsections corresponding to the IrReps of 
2?3d". These are listed in the first column. At the be- 
ginning of each subsection the second column gives the 
order parameter (OP) for the IrRep in terms of the no- 
tation A = AabcTa^T^^ cTc and arranges these into a 
vector. Next to this are the generators of symmetries 
under which the order parameter is invariant (generators 
and A commute). This completely characterizes one di- 
mensional represntations. 

In the case of the multi-component representations, 
particular values of the order parameter may have higher 
symmetries than the generic values - these are the dis- 
tinct phases. The values of the order parameter that 
produce the phase are given according to the vector rep- 
resentation in the second column. Next to these in the 
third column are the additional residual symmetries un- 
der which the phase is invariant, and the phase is labeled 
with the additional subscript. 

For example, let us take in section G in the first ta- 
ble. The second column of first line defines a vector 
for the representation. The third column states that 
all vectors in that representation are invariant under 
the time T time reversal operation. The next line says 
that when the vector is proportional to (1,0,0,1), i.e. 
A ex T^^T^^ '^'''v^'^v^ ' ^^^ symmetry is higher. The 
higher symmetries are in the third column; in particular, 
this vector is invariant under the combined G -\-t rotation 
and the Rc ■ Rt refiection in addition to the T rotation. 
The next line of the table says that when the vector takes 
the value (0, 0, 0, 1) the state is invariant under the refiec- 
tions Rt,Rc- Note that each phase is generally defined 
by a coset of values of the order parameter, which are 
all invariant under conjugate groups. We only list one 
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representative from each coset. For example, in the case 
of G, the (1, 0, 0, 1) vector is part of the coset of vectors 
(cos0,sin6', — sin6',±cos6'), 9€ [0,27r]. These are all in- 
variant under subgroups conjugate to the one listed in 
the third column. 

We use a, /3 as arbitrary real parameters when there 
is a continuous manifold of cosets. When listed under 
symmetries the symbols t, C , g and Sz mean the phase 
is invariant under the entire U{1) group generated. The 
symbol S means the phase is invariant under all spin ro- 
tations. There are several symmetry operations involv- 
ing rotations by tt or tt/2 in one the C/(l) groups. Like 
the spin reflections, these do not distinguish any of the 
phases so we do not list them. Five of the phases be- 
longing to the G representation are illustrated in Fig. [16] 
Notice that, according to the Landau theory, the tran- 
sition to the G -type phase can not occur directly but 
rather through the pattern with incommensurate period- 
icity. 




1. Normal Phases 

The normal phases are by definition invariant under 
spin and gauge transformation so we will suppress them. 
A product of Pauli matrices acting in different sub-spaces 
should be understood as a direct product. 



FIG. 16. (Color onlme)Sketch of the symmetric phases mag- 
netic and normal phases transforming according the G repre- 
sentation. (1) Gi normal state; (II) G2 normal state; (III) G3 
spin state; (IV) G4 spin state; (V) G5 spin state. Because of 
the absence of spin-orbit coupling the overall direction of the 
spins is arbitrary. 
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2. Magnetic phases 

We restore the spin symmetries but continue to sup- 
press the gauge symmetry. The high symmetry axis is 
arbitrarily chosen to be the z direction. 
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3. Superconducting Phases 

The order parameter M is defined by the non-zero ex- 
pectation values {ip''^ MTip'''). This M is listed under OP. 
M must contain an even number of Pauli matrices to be- 
cause of fermion anticonimutivity but may take complex 
values. 
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